Complex phases in the doped two-species bosonic Hubbard Model 
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We study a two-dimensional bosonic Hubbard model with two hard-core species away from half 
filling using Quantum Monte Carlo simulations. The model includes a repulsive interspecies inter- 
action and different nearest-neighbor hopping terms for the two species. By varying the filling we 
find a total of five distinct phases, including a normal liquid phase at higher temperature, and four 
different phases at lower temperature. We find an anti-ferromagnetically ordered Mott insulator 
and a region of coexistent anti-ferromagnetic and superfiuid phases near half filling. Further away 
from half filling the phase diagram displays a superfiuid phase and a novel phase inside the super- 
fiuid region at even lower temperatures. In this novel phase separated region, the heavy species 
has a Mott behavior with integer filling, while the lighter species shows phase separated Mott and 
superfiuid behaviors. 
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A major undertake in condensed matter physics is the 
study of exotic quantum phases. Strongly correlated sys- 
tems display quantum criticality, coexistent and inhomo- 
geneous phases, and secondary ordered phases close to 
quantum critical pointsii These exotic phenomena occur 
due to competing complex interactions of spin, charge, 
lattice, and orbital character. Competing phases in 
strongly correlated electronic systems^ occur in ruthen- 
ateSf^ manganiteS)^ cuprates;^ cobalt oxides?^ organic 
saltsi and heavy fermion materialsi^ Since competing 
phases can be tuned via external parameters, there are 
exciting opportunities for advance device applications of 
complex materials. 

Cold atoms experiments^ have become a playground 
for realizations of the HubbardiSiii and other strongly 
correlated model Hamiltonians, since model parame- 
ters can be tuned using laser and magnetic fields F^^ii^ 
Recently, there is an increasing interest in studies of 
mixtures of atoms due to the complexity associated 
with multiple species and the possibility of discovering 
novel phases. These include mixtures of bosonic and 
fermionic atomsi^ and mixtures of two different bosonic 
species.— Moreover, the experimental study of the 
^^K bosonic mixture in a three dimensional optical lat- 
ticei^ has motivated theoretical studies of the half-filled 
Hubbard model with two species with different masses or 
hopping integrals jiiiS These studies reveal a rich phase 
diagram. The concentration of carriers may also be con- 
trolled in experiments, which motivated our study. 

We consider two species of hard-core bosons a and b 
confined to a two-dimensional lattice. The Hamiltonian 
takes the form: 



where a\ {b\) and (bi) are the creation and annihilation 
operators, respectively, of hard-core bosons a (b), and 
n,° — ajcj , fi\ = b\b^ are the number operators of a, 
& bosons on site i, respectively. The sum runs 
over all distinct pairs of first neighboring sites i and j, 
ta(th) is the hopping integral between i and j sites for 
species a (6), and U"-^ is the strength of the interspecies 
repulsion. Since we are in the hard-core limit, creation 
and annihilation operators satisfy commutation rules on 
different sites i ^ j, [a, , a J = [a|,aj] ~ [a^ , aj] = 0, 
and anti-commutation rules on identical sites, {a|, a|} = 
{ a j , a j } =0 and ia^.aj] = 1 . 



We perform Quantum Monte Carlo simulations using 
the Stochastic Green Function algorithmic with Global 
Space-Time Updates^S for the canonical ensemble on 
L X L lattices with an inverse temperature /3 = 8L. 
Our results at half-filling reproduce the phase diagram 
of Ref. [l^. We focus on the unpolarized phase diagram, 
so our total density is p = N/L with N = Na+Nb = 2Na, 
with Na and TVf, the number of a and b particles, respec- 
tively. Since our intention is to explore the doping de- 
pendence of the model we restrict our simulation to the 
following parameters: ta ~ 0.08 it = t, and C/°^ = 6t, 
where t = 1. 

In Fig. [1] we look for signatures of ordering. To look 
for phase separation and Mott states, we calculate the 
chemical potential by adding one a and one b particle 
to the system as ^ = {E{N + 2) - E{N))/2. The su- 
perfiuid (SF) phase is detected by measuring the super- 
fluid density, psF, using the fluctuations of the wind- 
ing number, W , via Pollock and Ceperley's formula^i. 
The antiferromagnetic (AF) phase is characterized by 
a finite density-density static structure factor: S{k) ~ 

^ E] S2;p[ifc.(j5 — ^)] (n^^'^'nl"'^'). We find an AF phase 
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at half filling pi = 1. It is characterized by a vanishing 
compressibility, n = dp/d^ (top panel), a finite static 
staggered structure factor (bottom panel) , as well as AF 
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FIG. 1: (Color online) Top panel: The average density, p = 
N/L as a function of the chemical potential, fi. The vertical 
blue dashed line shows the Maxwell construction between pi 
and p2- The brown dotted-dashed lines show different phase 
boundaries at pa and p4 as discussed in the text. The inset 
shows a snapshot of the density profile at half filling, p = pi = 



I, with the blue (red) squares indicating 
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Bottom panel: The superfluid density (blue circles) and the 
staggered structure factor (green squares) as a function of 
p. The dot-dashed lines show different phase boundaries. All 
data are for L = 10, ^ = 80, ta = 0.08, U = 1.00 and U"'' = 6. 
Error bars are smaller than symbol sizes. 



ordering as shown in a snapshot of the density profile 
(inset of the top panel). Near half filling, p = N/L vs. 
fi displays a first-order phase transition between pi and 
P2 1.16 with instability, characterized by a negative 
slope, in a portion of the region. These two phases with 
densities pi and p2 coexist for any density value between 
the two end points. Since, in this region the system dis- 
plays finite values of S'(7r,7r) and psF we conclude that 
AF and SF phases coexist for any pi < p < P2- A ho- 
mogeneous SF state exists between p2 and p3 ^ 1.25 
identified by measuring the superfluid density psF- At 
P3 the superfluid density displays a decrease and the p 
versus p. plot shows a small bump. Another small feature 
is displayed in the top panel at p^ 1.52. Finally, the 
homogeneous SF phase continues until full filling. Next, 
we investigate the unexpected lowering of the superfluid 
density between p^ and p^. 

Fig. [2] shows snapshots of the average local density of 
both species for p^ < p = 1.44 < ^4. Simulations with 
open and periodic boundary conditions show clear evi- 
dence of phase separation. Species a always shows Mott 




FIG. 2: (Color online) Snapshot of the average local density 
of a and b particles for L = 20 and p = 1.44. Top panel: Open 
boundary conditions. The a and h local densities, (n°) and 
(n^), are shown on the left and right panels, respectively. For 
a particles, sites close to the boundary (red) have (nl) ~ 0, 
while the occupation of the central (blue) region is (n") ~ 1. 
For h particles (right panel) the density close to the boundary 
(blue region) is {n\) ~ 1, and at the center (green region), 
{n'l) ~ 0.60. Middle panel: The same quantities shown in the 
top panel but with periodic boundary conditions. Bottom 
panel: At the left the homogeneous density distribution of 
both a and b particles for p = 1.72. At the right, a sketch 
of the density profile for a simulation with periodic boundary 
conditions. 



behavior with ^ or 1, while species h displays 

Mott {n\) ~ 1 and SF coexistent phases. We can under- 
stand the tendency of the system to form such a mixed 
phase by extending the bosonic mean-field formalisroiSi^ 
to two-specics bosons in the hardcore limiti^i Using the 
Gutzwiller variational approach)^ the most general site 
factorized wave function can be written as 



*-n*^ = n [sill ^(sinfal+ cos 

i i 

+ cos^(sin^ + cos^al6l)]|0), (2) 



where 9, a and /? are variational parameters and |0) is 
the vacuum state. The energy function per lattice site 
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takes the form 



E 
Z2 



-ta sm cos y — - — j ~ tb sm sm" ' 



/3 



-C/^^cos^ -cos^ 
2 2 



(3) 



We minimize E subject to constraints on the densities 



r2 2^ 



and nb = j^ 



L2 ^2 ° ^2 ^ 2.2 ■ 

We solve these nonlinear equations by constructing an er- 
ror function and minimizing it. For p = 1.44, within this 
approximation, we find that the phase separated state is 
lower in energy than an homogeneous supcrfluid phase. 
The total energy of a supcrfluid with homogeneous den- 
sities pa = pb = 0.72 is 804.33 for L = 20. The bottom 
right panel of Fig. [5] represents the phase separated state 
as a central cross with area (2/i -I- Z)^ — Ah^ and four cor- 
ner squares of area h? . If we assume that in the cross 
= 1.0 and n\ = 0.60 for all i sites the mean-field en- 
ergy is 754.19. At the four squares n° = 0.0 and = 1.0 
and the total energy is zero since it is a Mott state. We 
can hypothesis that this new PS region becomes stable 
to reduce the potential energy. Consequently there is a 
critical U^^ above which the phase separated state is sta- 
ble. We estimate Uf - 4.8 for our chosen parameters, 
ta = 0.08, tb = 1.0 and p = 1.44. 
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FIG. 3: (Color online) Top panel: Average density, p, ver- 
sus chemical potential, /i, for different temperatures. Bottom 
panel: The superfluid density, psF, as a function of p for dif- 
ferent temperatures. All data are for system size L — 10. 
Error bars are smaller than symbol sizes. 

Next we study the temperature dependence of p and 
psF- The top panel of Fig.[3]shows p versus p for a system 



of size L = 10 and different inverse temperatures. Since 
there is a clear signature of phase separation for /? = 6.5 
but not for /3 ~ 1.5, we can conclude that the critical 
AF temperature occurs between these two temperatures. 
Similarly, from the psp vs. p curves for different temper- 
atures (see bottom panel of Fig. [3]) , wc can conclude that 
for this cluster size the phase separated region between 
P3 and p4 appears for temperatures between /3 = 6.5 and 
13—12. We infer the phase diagram by appropriate scal- 
ing of our finite-size results. 
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FIG. 4: (Color online) Bottom panel: The temperature, T, 
versus doping, S, phase diagram with equal population of each 
species. At half-filling, the system is antiferromagnetic (AF) 
below T^^ ~ 0.116 and normal liquid (NL) at larger temper- 
atures. By increasing the doping a discontinuous transition 
from AF to superfluid (SF) phase occurs, and an AF/SF phase 
separated region develops between 5i = and S2 — 0.06. Be- 
tween Ss = 0.13 and ^4 = 0.28 phase separation (PS) occurs 
with the heavy species becoming Mott while the the lighter 
one displays regions with either Mott or superfluid behaviors. 
Top left panel: Scaling behavior of the static staggered struc- 
ture factor for the continuous transition from AF to NL at 
half filling, 5i = (corresponds to the filled blue diamond in 
the bottom panel). Inset shows the collapse near to critical 
temperature with the Ising-like critical exponents. Top right 
panel: Superfluid density as a function of temperature for 
difi'erent system sizes at S — 0.07 (filled red diamond). Inset 
shows the finite size scaling to find the SF critical tempera- 
ture in the thermodynamic limit for the continuous transition 
at S — 0.07. Data points are based on simulation results, lines 
are guides to the eye. 



Fig. m displays the temperature, T, vs. doping, S ~ 
— ° ^ — - — — , phase diagram. In the thermodynamic 

limit, the AF phase only exists at half filling = 
and very low temperatures. The top left panel shows 
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the scaling of the AF to normal liquid (NL) continu- 
ous phase transition at Si = 0. This transition belongs 
to the two-dimensional Ising universality class for which 
the static staggered structure factor scales as ^(Tr, tt) = 
- r/^)^!/''), where / is a universal scal- 
ing function, and /? and v are the critical exponents for 
the order parameter and the correlation length, respec- 
tively. The factor 2f3/i' — 1/4 in two-dimensional sys- 
tems. Therefore we can read the critical temperature at 
the point where the ^(Tr, 7r)L^/'' vs. T curves for differ- 
ent system sizes cross. For our parameters T^^ = 0.116, 
and S{n,TT)L^/'^ vs. (T - T^^)L curves collapse. The 
AF phase is represented by a red line ending on a blue 
diamond in Fig. S) As we illustrate in previous figures, 
near half filling, wc find a discontinuous transition from 
AF to SF phase and a phase separation region for dop- 
ing 6i — 0.0 ^ 5 ^ 62 = 0.06 (dark blue region in 
Fig. 13]). The boundary of the AF/SF phase separated 
region is found by a Maxwell's construction of the p vs. 
fi plots. The PS region inside the SF phase exists for 
S3 = 0.13 ^ S S4 = 0.28. For a given temperature we 
determine its boundaries by estimating the filling where 
PSF starts decreasing (pa in Fig. [1]) or stop increasing 
(p4). For the rest of dopings we encounter a SF phase 
at low temperatures and an NL at higher temperatures. 
The top right panel of Fig. 2] shows the superfluid density 
as a function of temperature for different system sizes. 
The order parameter, superfluid density, has the univer- 
sal jump of = — at the critical pointi^ The tran- 

T TT 

sition from SF to NL is continuous and belongs to the 
Kosterlitz Thouless universality class. We find Tc in the 
thermodynamic limit by using the relation between the 
size dependent critical temperature Tc{L) and the cluster 

size: TJL) - TJoo) cx — ——^ The inset on the right 

top panel displays this scaling. For S — 0.07 {p = 1.14) 
wc find Tc ~ 0.254. Scaled transition points are shown as 
red diamonds in Fig. U) There is one triple point where 
AF, SF and NL phases coexist. 



In summary, by changing doping we find complex 
phases in the two-dimensional two-species hardcore 
bosonic Hubbard model for equal populations and un- 
equal hopping. Wc find a first order phase transition 
between the AF phase at half filling and a SF phase near 
half filling with a region of SF and AF coexistence. For 
a broad region of fillings and temperatures away from 
half filling, a SF phase is found. However, within the 
SF phase at finite doping we find a dome-shaped region 
containing an inhomogeneous phase. Density profiles of 
this novel phase separated region show that the heav- 
ier species displays Mott insulating behavior, while the 
lighter species is phase separated in Mott and superfluid 
regions. Farther from half filling, both species display an 
homogeneous superfluid phase. This complex phase di- 
agram reminds us about the phase diagram of cuprates. 
In particular, our phase diagram displays a region that 
is similar to the so-called "superconducting dome" away 
from half-filling. Further, we believe our work will en- 
courage experimental studies of this model on cold atoms 
traps. Indeed the experimental realization of half-filling 
AF phase is hard due to the low entropy associated with 
this phase, while the complex mixed phase that we iden- 
tify away from half-filling can be expected to have a high 
entropy and, then, be easier to obtain. To further ex- 
plore this complex phase diagram we are planning to ex- 
tend our simulations to polarized systems with different 
species population. 
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